Multiresolution Elastic Medical Image Registration in Standard Intensity Scale 
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Abstract 

Medical image registration is a difficult problem. Not 
only a registration algorithm needs to capture both large 
and small scale image deformations, it also has to deal with 
global and local image intensity variations. In this paper 
we describe a new multiresolution elastic image registration 
method that challenges these difficulties in image registra- 
tion. To capture large and small scale image deformations, 
we use both global and local affine transformation algo- 
rithms. To address global and local image intensity vari- 
ations, we apply an image intensity standardization algo- 
rithm to correct image intensity variations. This transforms 
image intensities into a standard intensity scale, which al- 
lows highly accurate registration of medical images. 



1 Introduction 

Image registration is important in many imaging appli- 
cations. For example, in diagnostic imaging, there is often 
a need for comparing two images for disease diagnosis or 
longitudinal studies. There is also a frequent need for regis- 
tering two images of different imaging modalities, e.g., MR 
and PET images. 

Image registration involves the development of a reason- 
able transformation between a pair of images, the source 
and target, such that the similarity between the transformed 
source image (registered source) and target image is opti- 
mum. The similarity measure should capture both large and 
small scale deformations (also known as displacements), 
together with global and local variations of image intensi- 
ties. Based on the nature of the transformation, registration 
methods can be categorized as rigid, affine and elastic. In 
rigid registration, the transformation includes global rota- 
tion and/or translation parameters. For Affine registration 
in particular, scaling parameters are also included. Registra- 
tion is considered elastic(deformable) if the transformation 
is able to express both global and local deformations. For 
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surveys of image registration including nonlinear medical 
image registration, see l^[T3l[T2ll8l[T8l[Tl[T4l. 

Although rigid and affine transformations are able to 
align images, they can only handle global deformations. In 
rigid registration, the recovered transformation itself has no 
clinical significance, however, in nonrigid registration the 
recovered transformation may have clinical significance [ 8 ] . 
Since motion and deformation characteristics are necessary 
for quantification of changes between images, transforma- 
tion should be found as accurate as possible. Except for 
a few studies [15] [19] [161 [TOl, most of the elastic defor- 
mations based nonrigid registrations rely on the assump- 
tion that image intensities remain constant between im- 
ages J2j [3] [12] Q~8], which is not always true and affects 
the accuracy of motion and deformations obtained from the 
transformation. 

To address this problem, a locally affine but globally 
smooth transformation model has been developed in the 
presence of intensity variations in fl5l . In addition to 6 
and 12 affine parameters for 2D and 3D registrations respec- 
tively, two more affine parameters are used to capture inten- 
sity variations during registration. In order to remove inef- 
ficiency and inaccuracy arising from certain circumstances, 
such as low-resolution images, Bayesian based importance 
sampling technique with the same spatially varying param- 
eters are used in lfl9l . In lfl6l . voxel based similarity mea- 
sures, such as normalized mutual information, are com- 
bined with B-spline based nonrigid transformation called 
free-form deformation (FFD). Since the intensity and con- 
trast between the pre- and post-contrast enchanced images 
vary, voxel based similarity measures are used because it 
is insensitive to intensity changes as a result of contrast 
enchancement. However, there is a trade-off between ac- 
curacy and computation time of FFD-based method. The 
local flexibility and computational complexity of the local 
motion model is related to the resolution of B-spline con- 
trol points. More control points may improve the registra- 
tion accuracy, but the computation time will also increase 
dramatically (6). In iflOl . a fast elastic multidimensional 



intensity-based image registration with a parametric model 
of deformation is presented. Although adding landmarks 
controls the smoothness of deformation field and using a 
multiresolution approach for both the image and the defor- 
mation model makes registration algorithm robust and fast, 
global solution of the optimization function cannot be guar- 
anteed due to manual identification of landmarks. 

In this paper, we present a multiresolution elastic image 
registration framework on images in the standard intensity 
scale. The standard intensity scale is obtained by a stan- 
dardisation procedure which corrects image intensity varia- 
tions |4|. In the standard scale, similar intensities will mean 
similar tissue properties. 

The paper is organised as follows: a detailed description 
of elastic registration used in this study is given in Section[2] 
A brief explanation of intensity standardization method is 
presented in Section [3] Multiresolution framework is ex- 
plained in SectionS] Experiments and promising results are 
given in Section|5]and concluding remarks in Section|6] 

2 Local Affine Transformation 

For 2D image registration, an affine transformation has 
six parameters, which can be determined if the coordinates 
of at least three non-colinear corresponding points in the 
images are known 0. As v = [x y 1] T represents the 
homogeneous spatial coordinates, let F(v, t) and F(v* , t — 
1) be the source and target images respectively. General 
affine transformation between source and target image can 
be modelled locally as: 
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F(v,t) 
1 



= F(v*,t-l) =F{A.v,t-l) 



(1) 

where t is the time, aj and as are parameters describing 
intensity variations between source and target images and 
A is affine transform matrix defined as: 



where D denotes a small spatial neighborhood. From a Tay- 
lor expansion of Eq ©, we obtain linear quadratic error 
function to be minimized: 
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where I is 3x3 identity matrix, / = [f x f y f t ] T and 

fx = fx(v,t) fy = fy(v,t) ft = ft(v,t) (5) 

are partial derivatives of image F on D. Open expression 
for the gradient based constraint equation can be ex- 
pressed further as: 

= ^2 (x.fx-ai + y-f x -a>2 + fx-a-3 + x.f y .a 4 

x,y£D 

+y-fyO,5 + fy-0.6 ~ f T -v) 2 

(6) 

where a = [a x a 2 a 3 04 a$ clq] t . More 
compact form can be obtained by defining f2 = 
[x-fx V-fx fx x.fy y.f y ,/ y ] T asin H3]: 



E(a) = ^2 lti T -a- f T -v] 
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Quadratic error function in Eq (0 can be minimized ana- 
lytically by differentiating it with respect to the unknown 
parameters a 



dE{a)/da= 2 ^ 
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Setting Eq (© equal to zero and solving for a parameters 
yields: 
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In general notation of affine transformation, 07 and ag con- 
trols the intensity variations between image pairs. In pro- 
posed registration algorithm since possible variations in in- 
tensities are captured by standardization method, Eq (HJ is 
minimized by setting a-j = 1 and ag = respectively. A 
simple way to estimate affine matrix parameters is to min- 
imize quadratic error function E(A) which can be defined 
as: 
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Since the velocity field at each image point has two compo- 
nents while the changes in image brightness at a point in the 
image plane due to motion yields only one constraint, the 
optical flow cannot be computed at a point in the image in- 
dependently of neighboring points without introducing ad- 
ditional constraints [9|. This additional constraint is based 
on the smoothness of parameters over domain D such that 
neighboring points on the domain D have similar velocities 
and the velocity field of the brightness patterns in the image 
varies smoothly almost everywhere. One way to express 



this additional constraint is to minimize the magnitude of 
the gradient of the optical flow velocity, which is: 
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where magnitude of A,; reflects the influence of smoothness 
term. Hence, the problem is to minimize the sum of the 
errors in the equation (0 and (fTOb . To obtain local affine 
parameters Si, we can differentiate E to tai{a) = E{a) + 
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(a) and set the result to be zero B9J. Solution for 
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where A is 6x6 diagonal matrix with elements A^, and is 
the component-wise average of a over domain D. Starting 
with the initial guess (a ) Q , at each step the next local 
affine paramters are computed and resultant system of linear 
equations are solved accordingly. 

3 Standardization of MR Image Intensity 
Scale 

Since MR image intensities do not have a fixed tissue 
meaning in image scale even within the same protocol, for 
the same body region, for images obtained on the same 
scanner, for the same patient, there is a strong need to trans- 
form image scale into standard intensity scale in order so 
that for the same MR protocol and body region, similar in- 
tensities will have similar tissue meaning [4|. 

Standardization is based on mapping image intensity his- 
tograms into a standard histogram. The method consists of 
two steps called training and transformation. In the training 
step, a set of images of the same body region are given as 
input to "learn" histogram-specific parameters, called land- 
marks. In the transformation step, any given image is stan- 
dardized with estimated histogram-specific landmarks ob- 
tained from the training step. 

3.1 Methods 

Based on the study [ 1 1 1, image histograms of the same 
body region are always of the same type. There are differ- 
ent types of histograms, unimodal, bimodal and generaliza- 
tion of both. Since most of the protocols studied in ifTTl 
produce bimodal histograms, bimodal histogram distribu- 
tion is used to extract histogram specific landmarks. In bi- 
modal histograms, one of the histogram specific landmarks 
is mode(fi) representing main foreground object in the im- 
age, as depicted in Figure 1 . 



1 aP can be obtained through the same equation by putting A 
null, 1151 



Figure 1. Location of the histogram spe- 
cific landmarks, =minimum gray value, 
m 2 =maximum gray value 




Other histogram specific landmarks denoted by p\ and 
P2 are extracted according to range of intensity of interest 
(IOI). Landmarks pi and p2 are defined according to mini- 
mum and maximum percentile values, pc\ and pc2, that are 
used to select IOI. 

In the training step, for image j, the landmarks pij, 
P2j and /ij obtained from the histogram of each image are 
mapped to the standard scale by mapping intensities from 
[Pij > P2j] to [si , S2] where si and S2 are minimum and max- 
imum intensities on the standard scale respectively. The for- 
mula for mapping x € [pij, P2j] to x' is the following ATI . 



X = Si + 



Pij 



P2j ~ Pij 



(s 2 - Si) 



(12) 



Figure 2 shows two separate linear mappings, the first 
from [pujfii] to [si,/i s ] and the second from [fii,p2i] 
to [^ s ,s 2 ]. Overall mapping, Tj(x), from [to^, TO 2 j] to 
\s'm, Sn„l can be summarized as follows: 



T i\ x ) = \ f _ \ 

\H S + (x - (^z^Jl if^<x<m 2l 

(13) 

where [7] converts any number yg 3? into closest integer 
Y such that Y > y or < y. Further details can be found 
in DH. 

3.2 Choosing the Standardization Param- 
eters 

Based on the experiments in lfTTl l4l. minimum and maxi- 
mum percentile values are set to pc\ = and pc2 = 99.8 re- 
spectively. On the standard scale, si and s 2 are set to si = 1 



Figure 2. The intensity mapping function for 
the transformation step 



Figure 3. Histograms before and after stan- 
dardization 




and s 2 = 4095. Figure 3 shows the histograms before in- 
tensity mapping and after intensity mapping respectively. It 
is easily seen that histograms are more similar in shape, lo- 
cation and distribution after standardization than before. It 
means that intensities have tissue-specific meaning after the 
standardization. 

Images in the first row of Figure 4 shows a source and 
target images in image scale. In the second row on the same 
figure, the same slices are displayed after standardization 
using the parameters defined above. 
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4 Multiresolution 



Image registration can be highly nonlinear and therefore 
many iterations may be required to reach a solution. 
An important method to reduce the amount of com- 
putational cost and deal with nonlinearities is to use a 
multilevel(multiresolution) image pyramid. Multilevel 
continuation is well established for optimization problems 
and systems of non-linear equations Q. As common to 
many other nonrigid registration algorithms, the regis- 
tration method we use includes two steps. After global 
registration has been done in the first step, locally affine 
globally smooth elastic registration on standard scale is 
performed. In order to achieve low computational cost and 
accelerate the registration process, coarse-to-fine strategy 
is used. Global affine registration is first performed at 
the coarsest level where convergence is fast because there 
are few data. The initial condition at the coarsest scale is 
arbitrary. Moreover, it is likely that the criterion to optimize 
has a reduced number of local optima; this is due to a loss 
of image details and results in enhanced robustness ifTTl . 



As seen from Figure 5, registration in finer level is per- 
formed with the result of the previous level as initial condi- 
tion. This process continues until the finest level is reached. 
Since the number of iterations performed at the finest level 
is more relevant than the other levels, for the computational 
cost of the whole optimization, it is very important that the 
initial condition for this last level be the best possible in or- 
der to reduce the amount of refinement necessary to reach 
convergence. To get optimal starting conditions, it is cru- 
cial that the coarse levels of the pyramid most represent the 
finest level |fl7l . In addition to this, the use of interpolation 
models (except linear interpolation) will change the inten- 
sity histogram of image after each warping. Available inter- 
polation methods vary in their computational complexity, 
speed and accuracy. To ensure a more accurate solution, we 
perform standardization after each warping/interpolation. 
Either small or large, intensity changes caused by the in- 
terpolation are captured by standardization. 



Figure 4. images in non-standard(1st row) 
and standard scale(2nd row) 



Figure 5. Coarsest-to-Finest Image 
Representation-Gaussian Pyramid 




5 Experiments and Results 

The registration algorithm is tested in two different ex- 
periments. The first experiment is global affine registration 
and the second is elastic registration, both in the multireso- 
lution framework. Accuracy of registration results is eval- 
uated quantitatively and qualitatively. Quantitative evalua- 
tion includes mean square differences in intensity between 
the pair of images |5 For visual assesment of registra- 
tion, a checkerboard image is created by picking alternating 
squares from the image pair, the registered source and tar- 
get. Spatial alignment can be investigated visually by look- 
ing at continuity of tissue boundaries in the two interleaved 
images. The alignment is good if there is no discontinu- 
ity of contours and the tissues merge smoothly across the 
checkerboard borders. 

5.1 Experiment I 

In the first experiment global deformations are captured 
with affine transformations in the multiresolution frame- 
work. For both the source and target images, four level 
coarsest to finest Gaussian pyramid is used. In each mul- 
tiresolution level, an optimal solution is determined and 
used as the starting point for the next level as seen in 
Figure 5. In each level, a single global affine transform 
is estimated with domain D, which is defined to be en- 
tire image. Transformed images are interpolated with cu- 
bic splines. Since intensities of images are changed in 
each warping/interpolation pair, image intensities are re- 
standardized into standard scale from image scale. 

2 Mean Square Error 




The algorithm is performed both on image scale and on 
standard intensity scale. Thirty deformed brain images are 
generated by applying random deformations to the original 
target images. Original source images are tried to be regis- 
tered to each deformed target images on standard intensity 
scale and image scale. Random deformations are obtained 
by randomly choosing parameters for affine matrix A. The 
resulting deformation field is normalized so that r.m.s dis- 
placement is at most 12 pixels. An example for affinely 
warped images is shown in Figure 6. 



Figure 6. Images are randomly deformed by 
affine transformations 




In Figure 7, three example registration results of ran- 
domly and affinely warped images are shown. The result- 



ing images clearly show that registered source images are 
in good agreement with target images. Registration quality 
is measured over 30 randomly deformed images by mean of 
the square of the differences in intensity (MSE). Experiment 
has been done both in image scale and on standard scale to 
show improvement in MSE sense. Table- 1 shows the MSE, 
maximum MSE and minimum MSE over 30 registration ex- 
amples on image scale and on intensity scale respectively. 



Figure 7. Resulting registration of images 
with random affine warps. Each row includes 
source, target and registered source 
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Table 1. Affine registration evaluations 



Scale/Error 


MSE 


Max MSE 


Min MSE 


On Image Scale 


0.1211 


0.1703 


0.0993 


On Standard Scale 


0.1060 


0.1457 


0.0857 



5.2 Experiment II 

In the second experiment, in addition to global de- 
formations, local deformations are captured with locally 
affine transformations in the multiresolution framework. 
The experiment starts with global affine registration in 
which pre-alignment of the source and target images is 
obtained. Similar to the first experiment, standardization 
is performed to handle image intensity variations in every 
warping/interpolation pair. The resulting global affine pa- 
rameters are used in the coarsest level of the Gaussian im- 
age pyramid to estimate local affine parameters on domain 
D, where D is now 5x5 pixels. Estimated local affine param- 
eters are used to transform source image in the next level of 



the Gaussian image pyramid. Accumulating all the trans- 
formations on each Gaussian image pyramid level yields a 
single final transformation that is able to capture both large 
and small scale motions. 

As in the previous experiment, the algorithm is per- 
formed both on image scale and on standard intensity 
scale. Thirty deformed brain images are generated by ap- 
plying random nonlinear deformations to the target images. 
Random nonlinear deformations are obtained according to 
equation given below: 

x 1 = m.x+ (-l) n2 .e na .sin(y/n 4 ) 

y' = n 5 .x + (-l) n °.e n ?.cos(y/n 8 ) (14) 

where [x, y], \x\ y'} are initial and transformed spatial coor- 
dinates of the images respectively, and, m, ..n$ are chosen 
randomly such that r.m.s displacement is at most 12 pix- 
els. An example for nonlinearly warped images is shown in 
Figure 8. 

Figure 8. Images are randomly deformed by 
nonlinear transformations 




In Figure 9, three example registration results of ran- 
domly and nonlinearly warped images are shown. Captur- 
ing signal intensity variations during registration process 
with intensity standardization method leads to assesment 
of visual comparision of registered source and target im- 
ages with warping grid. Evaluation of the registration re- 
sults is summarized in Table-2. The table shows that large 
and small scale deformations are captured accurately on the 
standard intensity scale. Resulting images have fixed in- 
tensity meanings even there is large intensity variations ini- 
tially. 

The resulting registered images and deformation fields 
show that standardization of intensity scales improves the 
accuracy of registration. 



Figure 9. Resulting registration of images 
with random nonlinear warps. Each row in- 
cludes source, target, registered source and 
estimated warping grid 




Table 2. Elastic registration evaluations 



Scale/Error 


MSE 


Max MSE 


Min MSE 


On Image Scale 


0.0204 


0.0356 


0.0155 


On Standard Scale 


0.0194 


0.0320 


0.0148 



Another method to evaluate proposed registration 
method is visual examination of checkerboard images. Fig- 
ure 10 shows an elastic registration example together with 
checkerboard image illustrating how well the image pair is 
registered. Checkerboard image includes white and black 
squares corresponding to intensity values taken from the 
registered source and the target image respectively. Our 
overall observation from experimental results is that mul- 
tiresolution elastic registration on standard intensity scale 
can capture both local and global deformations with high 
accuracy. 

6 Conclusion 

In this paper, a multiresolution elastic medical image 
registration is presented. The multiresolution framework 
leads to a robust and fast algorithm. Variations in image 
intensity histograms at each pyramid level is corrected by 
the intensity standardization method in order to provide an 
efficient framework for registration. The standardization 
method enables similar image intensities mean similar tis- 
sue contents, which leads to less parameters in registration. 
Qualitative and quantitative evaluation of experimental re- 
sults indicate that small and large scale deformations are 



Figure 10. First row includes the source and 
target images, second row shows registered 
source with warped grid and checkerboard 
image for visual assesment 




captured with high accuracy. The examples presented here 
use 2D images but the proposed algorithm is valid in 3D as 
well. 
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